! **************************************************************************************
!  GridHDF5 (version 2.1) 
!  Rob Crain (2007), Institute for Computational Cosmology, Durham
!  Code produces hdf5 grids of smoothed quantities from Gadget-3 snapshots
! **************************************************************************************

program main
  use logs
  use my_mpi
  use hdf5_wrapper
  use arg_module, only : isnap,iframe
  use camera_module
  use parameter_module
  use transform_module
  use myfiles
  use smoothing_variables, only : nvar_by_type,var
  use grid
  use particles
  use header
  use units_and_constants
  implicit none
  integer*4 :: itype,ivar, ierr
  character(len=16) :: ctype,cPE
  real*8, external :: compute_comoving_distance  
  

  ! Perform startup tasks and read the input files
  call startup()
  call parse_paramfile()
 
  ! Are we using more cores than z-slices? If so, terminate
  if(NumPEs > params%nz) call abort("NumPEs must be smaller than or equal to nz")

  ! Are we using a camera file?
#ifdef USE_CAM_FILE
  call read_camerafile()
#endif

  write(logunit,*)' -MyPE = ',MyPE
  call flush(logunit)

  ! Dump some useful info to the screen
  write(logunit,*)' -camera details'
  write(logunit,*)'  -isnap  =',isnap
  write(logunit,*)'  -L      = (',L(1),'|',L(2),'|',L(3),')'
#ifdef CUBE
  write(logunit,*)'  -n      = (',params%nx,'|',params%ny,'|',params%nz,')'
#else
  write(logunit,*)'  -n      = (',params%nx,'|',params%nz,')'
#endif
#ifdef CUBE
  write(logunit,*)'  -L/pix  = (',L(1)/real(params%nx),'|',L(2)/real(params%ny),'|',L(3)/real(params%nz),')'
#else
  write(logunit,*)'  -L/pix  = (',L(1)/real(params%nx),'|',L(3)/real(params%nz),')'
#endif
  write(logunit,*)'  -Centre = (',x0,',',y0,',',z0,')'
#ifdef VELOCITY
  write(logunit,*)'  -V_sys  = (',vx0,',',vy0,',',vz0,')'
#endif
#ifdef USE_EIGENVECTORS  
  write(logunit,*)'  -xprime = (',real(xprime(1)),'|',real(xprime(2)),'|',real(xprime(3)),')'
  write(logunit,*)'  -yprime = (',real(yprime(1)),'|',real(yprime(2)),'|',real(yprime(3)),')'
  write(logunit,*)'  -zprime = (',real(zprime(1)),'|',real(zprime(2)),'|',real(zprime(3)),')'
  call flush(logunit)
#else
  write(logunit,*)'  -Angle (deg)   = (',theta,',',phi,',',psi,')'
  call flush(logunit)
#endif  

  ! Main input routine
  call read_variables()
  
  ! Separate smoothing processes by particle types
  do itype=0,5
     
     ! Report which species we're playing with
     if(nvar_by_type(itype) > 0)then 
        write(ctype,'(I1.1)')itype
        write(logunit,*)' -Operating on PartType'//trim(adjustl(ctype))
        do ivar=1,nvar_by_type(itype)
           write(logunit,*)'  -var = ',trim(var(itype,ivar))
        enddo
       call flush(logunit)

#ifdef FORCE_EXTERNAL_HSML
       ! In this case we'll always look for smoothing lengths in an external file
       call read_data(itype,.false.) 
       call read_hsmlfile(itype)
#else
       ! In this case, we'll use snapshot hsml if available
        selectcase(itype)
        case(0)
           call read_data(itype,.true.) 
        case default
           call read_data(itype,.false.) 
	   call read_hsmlfile(itype)
        endselect
#endif

        ! Compute angular<->physical scale
        ! NB: if z=0, Dcom=0, so the conversion is meaningless
        ! so set a minimum comoving distance of half the boxsize (in physical units)
        Dcom = compute_comoving_distance(Redshift,Omega0,OmegaLambda,HubbleParam)
        Dcom = max(Dcom,(Boxsize/(2.0d0*HubbleParam))) ! Mpc (no h)
        arcsec_per_radian = (1.0d0/3.6d3)*(Constants%PI/1.8d2)
        comoving_angular_scale = 1.0d3*Dcom*dtan(arcsec_per_radian)
        physical_angular_scale = ExpansionFactor*comoving_angular_scale
        write(logunit,*)' -comoving distance scale (Mpc)  = ',real(Dcom)
        write(logunit,*)' -comoving angular_scale (kpc/") = ',real(comoving_angular_scale)
        write(logunit,*)' -physical angular_scale (kpc/") = ',real(physical_angular_scale)
        call flush(logunit)
        
        ! Perform coordinate transform
        call coordinate_transform()

        ! Set the number of planes per process
        call set_xyz_grids(isnap)
         
        ! Determine if, and where, particles are to be exported 
        call build_exchange_lists(itype)
 
        ! Execute the exchange
        call exchange_particles(itype)

        ! Exchange copies of ghost particles that contribute to external grids
        call exchange_ghosts(itype)
 
        ! Allocate the local smoothing grid
        call allocate_local_grids(itype)           

	! Smooth to the grid 
        call sph2grid(itype)

        ! Reduce the grids onto root process and write HDF5 file 
#ifndef NONEVOLVING
        iframe = isnap
#endif
        call mpi_barrier(mpi_comm_world,ierr)
        call write_hdf5_data(isnap,iframe,itype)
        
        ! Free up memory
        call deallocate_local_grids()

     endif
  
  enddo
     
  call flush(logunit)
  call mpi_barrier(mpi_comm_world,ierr)
     
  call endrun()

  stop
end program main

